home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_4.5 / XBDY / PROC_POL.C < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  11.8 KB  |  436 lines
  1. /* 
  2.  * proc_pol.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * PROC(ess)_POLY(gon)
  11.  *
  12.  * polygon analysis:
  13.  * -----------------
  14.  *   reconstruct polygon;
  15.  *   determine area by pixel counting
  16.  *            (and, in poly_moments(), by evaluating zero order moment);
  17.  *
  18.  *   evaluate curvature energy based on Zahn chain code
  19.  *
  20.  *   compute centroid, central moments and principal axes;
  21.  *
  22.  *   find 2**n new vertices equi-spaced along contour
  23.  *   determine radius vector relative to centroid; shift to zero mean;
  24.  *   reconstruct angular bend function with respect to new set of vertices;
  25.  *
  26.  *   compute power_spectrum and autocorrelation by AP FFT,
  27.  *           for angular bend function (-->tan mode) 
  28.  *           for radial function (-->rad mode)
  29.  *           for both (-->all)
  30.  *   save to file (writes file for plotting by routine acmp);
  31.  *
  32.  *   evaluate Zahn-Roskies Fourier descriptors
  33.  *   write them to file (writes file for plotting by modeplt, multiplt);
  34.  *
  35.  */
  36.  
  37. #include <stdio.h>
  38. #include <stdlib.h>
  39. #include <math.h>
  40. #include "ip.h"
  41.  
  42. #define    ON_MODE        1
  43. #define    OFF_MODE    0
  44. #define    C_HEIGHT    10
  45. #define    C_WIDTH        10
  46. #define    C_COLOR        255
  47.  
  48.  
  49. #define N2_LIMIT    64L
  50. #define    MAX_ZRFD    15L
  51. #define N_MOMENTS    15
  52. #define N_MODE_PARMS    2
  53.  
  54. #define ENLARGE        1.0
  55. #define NOLOOP         0
  56. #define MAX_ROW     2048
  57. #define LOOP_COMPLETE    8.0
  58.  
  59. #define SQ2        1.414213562
  60. #define ZERO        0.0
  61. #define FABS(a)        ((a) > ZERO ? (a) : - (a))
  62.  
  63. #define ON        1
  64. #define OFF        0
  65. #define SHIFT_CENTROID    ON
  66. #define    SAVE_SHAPE    OFF
  67.  
  68. #define    EVAL_CURV_EN    OFF        /* evaluate curv. energy */
  69. #define    POLY_FIT    ON             /* invoke Wall-Danielsson fit */
  70. #define    SPEC_ANAL    ON            /* perform spectral analysis */
  71. #define    ZAHN_ROSKIES    ON         /* evaluate Zahn-Roskies Fourier descript. */
  72. #define    DRAW_BDY_ONLY    OFF       /* draws reconstructed polygonal contour */
  73.  
  74. #define DISPLAY_DATA
  75.  
  76.  
  77. struct polygon *
  78. select_poly (struct polygon *poly_head, Image * imgIO, int value)
  79. {
  80.   struct polygon *current_poly;
  81.  
  82.   int c;
  83.   int retval;
  84.  
  85.   int rad_mode = OFF, tan_mode = OFF;
  86.  
  87.   float *s, *s_ap, *s_n2;
  88.   float *moments;
  89.   float *phi_k;
  90.   float *r_n2;
  91.  
  92.   float *n2_phi_star;
  93.   float *acf, *p_spec;
  94.   int *mode_parms;
  95.   float *a_n = NULL, *b_n = NULL;
  96.  
  97.   int i;
  98.   long nv, nv2;
  99.   long n_order;
  100.   int n_mom = N_MOMENTS;
  101.   int n_mp = N_MODE_PARMS;
  102.  
  103.   unsigned int pixel_count;
  104.  
  105.   double av_dirn = 0.0;
  106.   double c_len, ap_c_len, res_c_len;
  107.   double tot_phi;
  108.   double curv_en, e_bend, se_bend;
  109.   double p2a_ratio;
  110.  
  111.   struct Bdy bd, *bdp = &bd;
  112.   struct spoint *hullctr;
  113.   struct spoint cursor, *pc = &cursor;
  114.   struct spoint *v, *v_ap, *v_n2, *v_n2h;
  115.   struct spoint *pvc;
  116.   int hull_area;
  117.  
  118.   FILE *fOut, *fOut1;
  119.   char in_buf[IN_BUF_LEN];
  120.  
  121. /*
  122.  * cycle through list of boundary polygons
  123.  */
  124.   printf ("...enter n to cycle, y to select poly...\n");
  125.   current_poly = poly_head;
  126.   for (;;) {
  127.     pc->x = current_poly->first_x;
  128.     pc->y = current_poly->first_y;
  129.     printf ("poly starting at (%d, %d) select?(y/n) ", pc->x, pc->y);
  130.  
  131.     if ((c = readlin (in_buf)) == 'y')
  132.       break;
  133.     if ((current_poly = current_poly->next_poly) == NULL) {
  134.       printf ("At end of polygon list!!  Return to first polygon? (y/n)");
  135.       if ((c = readlin (in_buf)) == 'n')
  136.         return (current_poly);
  137.       else
  138.         current_poly = poly_head;
  139.     }
  140.   }
  141.  
  142.  
  143.   nv = current_poly->poly_points;
  144.   nv2 = 2L;
  145.   while (nv2 < nv)
  146.     nv2 *= 2L;
  147.   if (nv2 >= N2_LIMIT)
  148.     nv2 = N2_LIMIT;
  149.  
  150.  
  151. /*
  152.  * allocate memory 
  153.  */
  154.   if ((v = (struct spoint *) malloc ((nv + 1) * sizeof (struct spoint))) == NULL)
  155.       exitmess ("\n...mem allocation for v failed\n", 1);
  156.   if ((v_ap = (struct spoint *) malloc ((nv + 1) * sizeof (struct spoint))) == NULL)
  157.       exitmess ("\n...mem allocation for v_ap failed\n", 1);
  158.   if ((phi_k = (float *) malloc ((nv + 1) * sizeof (float))) == NULL)
  159.       exitmess ("\n...mem allocation for phi_k failed\n", 1);
  160.   if ((s = (float *) malloc ((nv + 1) * sizeof (float))) == NULL)
  161.       exitmess ("\n...mem allocation for s failed\n", 1);
  162.   if ((s_ap = (float *) malloc ((nv + 1) * sizeof (float))) == NULL)
  163.       exitmess ("\n...mem allocation for s_ap failed\n", 1);
  164.  
  165.  
  166.   if ((moments = (float *) malloc (n_mom * sizeof (float))) == NULL)
  167.       exitmess ("\n...mem allocation for moments failed\n", 1);
  168.  
  169.  
  170.   if ((v_n2 = (struct spoint *) malloc ((nv2 + 1) * sizeof (struct spoint))) == NULL)
  171.       exitmess ("\n...mem allocation for v_n2 failed\n", 1);
  172.   if ((v_n2h = (struct spoint *) malloc ((nv2 / 2) * sizeof (struct spoint))) == NULL)
  173.       exitmess ("\n...mem allocation for v_n2h failed\n", 1);
  174.   if ((r_n2 = (float *) malloc (nv2 * sizeof (float))) == NULL)
  175.       exitmess ("\n...mem allocation for r_n2 failed\n", 1);
  176.   if ((s_n2 = (float *) malloc ((nv2 + 1) * sizeof (float))) == NULL)
  177.       exitmess ("\n...mem allocation for s_n2 failed\n", 1);
  178.  
  179.   if ((n2_phi_star = (float *) malloc (nv2 * sizeof (float))) == NULL)
  180.       exitmess ("\n...mem allocation for n2_phi_star failed\n", 1);
  181.  
  182.   if ((acf = (float *) malloc (nv2 * sizeof (float))) == NULL)
  183.       exitmess ("\n...mem allocation for acf failed\n", 1);
  184.   if ((p_spec = (float *) malloc (nv2 * sizeof (float))) == NULL)
  185.       exitmess ("\n...mem allocation for p_spec failed\n", 1);
  186.   if ((mode_parms = (int *) malloc (n_mp * sizeof (int))) == NULL)
  187.       exitmess ("\n...mem allocation for mode_parms failed\n", 1);
  188.  
  189. /*
  190.  * evaluate arc_length and Zahn angular bend function
  191.  */
  192.   c_len = arc_length (current_poly->d_l_ptr, s, nv);
  193.  
  194.   tot_phi = angular_bend (current_poly->d_phi_ptr, phi_k, nv);
  195.  
  196. /*
  197.  * reconstruct polygon vertices
  198.  */
  199.   pixel_count = reconstruct_poly (current_poly->d_phi_ptr,
  200.                                   current_poly->d_l_ptr, nv, v, pc,
  201.                                   current_poly->first_edge_out,
  202.                                   current_poly->total_delta_phi,
  203.                                   imgIO, value);
  204.  
  205.   printf ("\n...object area (pixel count): %u", pixel_count);
  206.  
  207.  
  208. /*
  209.  * evalate curvature energy from chain code representation
  210.  * note:
  211.  *   to obtain robust estimates of curv_en, evaluation should really
  212.  *   not be attempted before smoothing, e. g. via polygonal approx.,
  213.  *   see below;
  214.  */
  215.   if (EVAL_CURV_EN == ON) {
  216.     curv_en = curvature_energy (current_poly->d_phi_ptr,
  217.                                 current_poly->d_l_ptr, nv);
  218.     printf ("\n...chain-code derived curv. energy: %f\n", curv_en);
  219.   }
  220.  
  221. /*
  222.  * Wall-Danielsson polygonal approximation
  223.  */
  224.   res_c_len = p2a_ratio = e_bend = -1.0;
  225.   if (POLY_FIT == ON) {
  226.     fill_bdy_structure (bdp, v, nv, tot_phi);
  227.     hullctr = poly_hull (bdp, v_ap, av_dirn, &hull_area, imgIO, GRAY);
  228.     printf ("\tlocation of convex hull center:");
  229.     printf (" (%3d,%3d)\n", hullctr->x, hullctr->y);
  230.  
  231.     if ((nv - bdp->an) > 0) {
  232.       printf ("\n...approx poly size: %ld\n", nv = bdp->an);
  233.  
  234.       v_ap = (struct spoint *) realloc (v_ap, (nv + 1) * sizeof (struct spoint));
  235.  
  236.       nv2 = 2L;
  237.       while (nv2 < nv)
  238.         nv2 *= 2L;
  239.     }
  240.  
  241. /*
  242.  * evaluate polygon moments
  243.  */
  244.     pvc = poly_moments (nv, v_ap, moments, imgIO, GRAY);
  245.     printf ("\n...centroid pos: (%3d, %3d)\n", pvc->x, pvc->y);
  246.  
  247. /*
  248.  * evaluate arc_length (partial sums of edge lengths) from set of
  249.  * vertices for approximated polygon;
  250.  */
  251.     ap_c_len = vert_to_clen (v_ap, s_ap, nv);
  252.  
  253. /*
  254.  * resample new polygon at 2**n sites
  255.  */
  256.     n2_vert (v_ap, s_ap, nv, v_n2, nv2, imgIO, 192);
  257.   }
  258.   else {
  259.     pvc = poly_moments (nv, v, moments, imgIO, GRAY);
  260.     n2_vert (v, s, nv, v_n2, nv2, imgIO, value);
  261.   }
  262.  
  263. /*
  264.  * evaluate contour length, global shape for resampled polygon;
  265.  */
  266.   printf ("\n...geom. parameters for polygon:\n");
  267.   res_c_len = vert_to_clen (v_n2, s_n2, nv2);
  268.   p2a_ratio = shape_parm (res_c_len, moments);
  269.  
  270. /*
  271.  * evaluate total bending energy, given by the contour integral of 
  272.  * the square of the curvature (see also spec_bend_en() in psaf.c)
  273.  * note:
  274.  *   oversampling (generating densely spaced vertices) amplifies
  275.  *   contour noise and leads to siginificant overestimates in
  276.  *   the bending energy (easily tested with circular shapes for
  277.  *   which the expected bending energy is given by 2*PI/R)
  278.  */
  279.   for (i = 0; i < nv2 / 2; i++)
  280.     *(v_n2h + i) = *(v_n2 + 2 * i);
  281.   e_bend = cont_bend_en (v_n2h, nv2 / 2);
  282.  
  283.   printf (" c_len: %lf, p2a_ratio: %lf, e_bend: %lf\n",
  284.           res_c_len, p2a_ratio, e_bend);
  285.  
  286. /*
  287.  * evaluate spectral densities
  288.  */
  289.   if (SPEC_ANAL == ON) {
  290.     printf ("\n...perform radial spectral analysis?(y/n)");
  291.     c = readlin (in_buf);
  292.     if (c == 'y')
  293.       rad_mode = ON;
  294.   }
  295.   *(mode_parms + 1) = SHIFT_CENTROID;
  296.   if (tan_mode == ON) {
  297.     *(mode_parms + 0) = tan_mode;
  298.  
  299.     vert_to_phi_star (v_n2, nv2, n2_phi_star);
  300.  
  301.     spec_and_corr (nv2, n2_phi_star, acf, p_spec);
  302.  
  303.     se_bend = spec_bend_en (nv2, p_spec, res_c_len);
  304.   }
  305.   if (rad_mode == ON) {
  306.     *(mode_parms + 0) = 1 - rad_mode;
  307.  
  308. /*
  309.  * evaluate radial shape function
  310.  */
  311.     r_dev (pvc, v_n2, r_n2, nv2, moments, imgIO, value);
  312.  
  313. #ifdef DISPLAY_DATA
  314.     printf ("\n...output of function r_dev():\n");
  315.     for (i = 0; i < nv; i++)
  316.       printf ("   r(%3d): %f\n", i, *(r_n2 + i));
  317. #endif
  318.  
  319. /*
  320.  * evaluate spectrum and (auto)correlation function of 
  321.  * radial displacement function r_n2
  322.  */
  323.     spec_and_corr (nv2, r_n2, acf, p_spec);
  324.  
  325.     se_bend = spec_bend_en (nv2, p_spec, res_c_len);
  326.   }
  327.  
  328. /*
  329.  * write spectral data to file (.acm)
  330.  */
  331.   c = 0;
  332.   printf ("\n...write file (.acm)?(y/n)");
  333.   while ((c != 'y') && (c != 'n'))
  334.     c = readlin (in_buf);
  335.  
  336.   if (c == 'y') {
  337.     printf ("...enter fn [drive:][\\directory\\]: ");
  338.     readlin (in_buf);
  339.  
  340.     if ((fOut = fopen (in_buf, "w")) == NULL) {
  341.       printf ("\n...cannot open output file %s", in_buf);
  342.       exit (1);
  343.     }
  344.  
  345.     write_acm_file (fOut, nv2, p_spec, acf, moments, n_mom,
  346.                     mode_parms, n_mp, c_len, (double) pixel_count,
  347.                     res_c_len, p2a_ratio, e_bend);
  348.  
  349.     fclose (fOut);
  350.   }
  351.   if (ZAHN_ROSKIES == ON) {
  352.     printf ("\n...compute Fourier descriptors(y/n)?");
  353.     if ((c = readlin (in_buf)) == 'y') {
  354.       for (;;) {
  355.         printf ("\n...enter # of Fourier descriptors: ");
  356.         readlin (in_buf);
  357.         retval = sscanf (in_buf, "%ld", &n_order);
  358.         if ((retval == 0) || (retval == EOF)) {
  359.           printf ("\n...try again: ");
  360.           printf ("...enter # of Fourier descriptors: ");
  361.         }
  362.         else if (n_order > MAX_ZRFD) {
  363.           printf ("\n...current max for FD: %ld\n",
  364.                   MAX_ZRFD);
  365.           printf ("\n...try again: ");
  366.           printf ("...enter # of Fourier descriptors: ");
  367.         }
  368.         else
  369.           break;
  370.       }
  371.  
  372. /* 
  373.  * allocate memory and evaluate Fourier descriptors
  374.  */
  375.       if ((a_n = (float *) calloc (MAX_ZRFD, sizeof (float))) == NULL)
  376.           exitmess ("\n...mem alloc for a_n failed\n", 1);
  377.       if ((b_n = (float *) calloc (MAX_ZRFD, sizeof (float))) == NULL)
  378.           exitmess ("\n...mem alloc for b_n failed\n", 1);
  379.  
  380.       msdescriptors (current_poly->d_phi_ptr,
  381.                      current_poly->d_l_ptr, nv, a_n, b_n, n_order);
  382.  
  383. /*
  384.  * write descriptors to file (.fd)
  385.  */
  386.       c = 0;
  387.       printf ("\n...write file (.fd)?(y/n)");
  388.       while ((c != 'y') && (c != 'n'))
  389.         c = readlin (in_buf);
  390.       if (c == 'y') {
  391.         printf ("...enter fn [drive:][\\directory\\]: ");
  392.         readlin (in_buf);
  393.  
  394.         if ((fOut1 = fopen (in_buf, "w")) == NULL) {
  395.           printf ("\n...cannot open output file %s", in_buf);
  396.           exit (1);
  397.         }
  398.  
  399.         write_fd_file (fOut1, n_order, a_n, c_len,
  400.                        (double) pixel_count, (float) p2a_ratio);
  401.  
  402.         fclose (fOut1);
  403.       }
  404.     }
  405.   }
  406.   printf ("\n-->polygon analysis completed<--\n");
  407.  
  408. /*
  409.  * deallocate memory 
  410.  */
  411.   free (v);
  412.   free (v_ap);
  413.   free (v_n2);
  414.   free (v_n2h);
  415.   free (s);
  416.   free (s_ap);
  417.   free (s_n2);
  418.   free (phi_k);
  419.   free (moments);
  420.   free (r_n2);
  421.  
  422.   free (n2_phi_star);
  423.  
  424.   free (acf);
  425.   free (p_spec);
  426.  
  427.   if (ZAHN_ROSKIES == ON) {
  428.     if (a_n != NULL)
  429.       free (a_n);
  430.     if (b_n != NULL)
  431.       free (b_n);
  432.   }
  433.  
  434.   return (current_poly);
  435. }
  436.